package jscinttools.utilities.scintutilities;

import hep.aida.*;
import jhplot.P1D;
import java.io.*;
import fr.in2p3.ipnl.aidautilities.*;
import java.lang.*;
//import java.lang.Object.*;
//import static java.lang.Math.*;


//import aidaCode.*;


/**
 *  Various Static Methods for Scintillation Analysis
 *
 *@author     Marc-Antoine Verdier
 *@created    6 décembre 2010
 */
public final class ScintUtil
{
	private final static IAnalysisFactory af = IAnalysisFactory.create();


	/**
	 *  Constructor for the ScintUtil object
	 */
	private ScintUtil() { }


	/**
	 *  Returns the light yield correction factor, based on the time window as well
	 *  as the decay times. Does not take rise into account (should not be a big
	 *  effect).
	 *
	 *@param  amp     The array with the amplitudes of the average event
	 *      components.
	 *@param  tau     The array of time constants
	 *@param  window  The time window (same units as time constants)
	 *@return         The correction factor (should be >= 1 by definition)
	 */
	public static double calcLYCorr(double[] amp, double[] tau, double window)
	{
		double num = 0;
		double denom = 0;
		for (int i = 0; i < amp.length; i++)
		{
			num += (amp[i] * tau[i]);
			denom += (amp[i] * tau[i] * (1 - Math.exp(-1.0 * window / tau[i])));
		}
		return num / denom;
	}


	/**
	 *  Returns the error on the light yield correction factor, based on the time
	 *  window as well as the decay times.
	 *
	 *@param  amp     The array with the amplitudes of the average event components
	 *@param  ampErr  The array of errors on the amplitudes of the average event
	 *      components
	 *@param  tau     The array of time constants
	 *@param  tauErr  The array of errors on the time constants
	 *@param  window  The time window (same units as time constants)
	 *@return         The error on the correction factor
	 */
	public static double calcLYCorrErr(double[] amp, double[] ampErr, double[] tau, double[] tauErr, double window)
	{
		double num = 0;
		double denom = 0;
		for (int i = 0; i < amp.length; i++)
		{
			num += (amp[i] * tau[i]);
			denom += (amp[i] * tau[i] * (1 - Math.exp(-1.0 * window / tau[i])));
		}
		double corrFacErr = 0;
		for (int i = 0; i < amp.length; i++)
			corrFacErr += Math.pow(((1.0 / num) - ((1.0 - Math.exp(-1.0 * window / tau[i])) / denom)) * tau[i] * ampErr[i], 2) + Math.pow(((1.0 / num) - ((1.0 - (1.0 + (window / tau[i]) * Math.exp(-1.0 * window / tau[i]))) / denom)) * amp[i] * tauErr[i], 2);

		corrFacErr = (num / denom) * Math.sqrt(corrFacErr);
		return corrFacErr;
	}


	/**
	 *  Description of the Method
	 *
	 *@param  compIndex  The index of the component
	 *@param  amp        The array with the amplitudes of the average event
	 *@param  tau        The array of time constants
	 *@return            Description of the Return Value
	 */
	public static double calcLYFrac(int compIndex, double[] amp, double[] tau)
	{
		double denom = 0;
		double num = amp[compIndex] * tau[compIndex];
		for (int i = 0; i < amp.length; i++)
			denom += amp[i] * tau[i];

		return num / denom;
	}


	/**
	 *  Description of the Method
	 *
	 *@param  compIndex  The index of the component
	 *@param  amp        The array with the amplitudes of the average event
	 *@param  ampErr     The array of errors on the amplitudes of the average event
	 *@param  tau        The array of time constants
	 *@param  tauErr     The array of errors on the time constants
	 *@return            Description of the Return Value
	 */
	public static double calcLYFracErr(int compIndex, double[] amp, double[] ampErr, double[] tau, double[] tauErr)
	{
		double theSum = 0;
		double theComp = amp[compIndex] * tau[compIndex];
		for (int i = 0; i < amp.length; i++)
			theSum += amp[i] * tau[i];

		double lyFracErr = 0;
		double sumErr = 0;
		for (int i = 0; i < amp.length - 1; i++)
			if (i != compIndex)
				sumErr = Math.pow(tau[i] * ampErr[i], 2) + Math.pow(amp[i] * tauErr[i], 2);

		lyFracErr = (theComp / theSum) * Math.sqrt(Math.pow(1.0 / theComp + 1.0 / theSum, 2) * (Math.pow(tau[compIndex] * ampErr[compIndex], 2) + Math.pow(amp[compIndex] * tauErr[compIndex], 2)) + Math.pow(1.0 / theSum, 2) * sumErr);
		return lyFracErr;
	}


	/**
	 *  Gets the index of the best chi_2 of an array of IFitResult
	 *
	 *@param  theResult  The array of IFitResult
	 *@return            The index of the best chi_2
	 */
	public static int getBestChi2Index(IFitResult[] theResult)
	{
		int index = 0;
		double chi2_0 = theResult[0].quality();
		for (int i = 0; i < theResult.length - 1; i++)
			if ((theResult[i + 1].quality() < chi2_0) && (theResult[i + 1].errors()[0] != 0.0))
			{
				chi2_0 = theResult[i + 1].quality();
				index++;
			}

		return index;
	}


	/**
	 *  Normalize the LY to 1 at RoomT
	 *
	 *@param  theDPS      the LY DPS
	 *@param  roomTLY     The Value of the RoomT LY
	 *@param  roomTLYErr  The Value of the error of the RoomT LY
	 *@return             The normalized DPS
	 */
	public static IDataPointSet normalizeToRoomT(IDataPointSet theDPS, double roomTLY, double roomTLYErr)
	{
		for (int i = 0; i < theDPS.size(); i++)
		{
			theDPS.point(i).coordinate(1).setValue(theDPS.point(i).coordinate(1).value() / roomTLY);
			theDPS.point(i).coordinate(1).setErrorMinus((theDPS.point(i).coordinate(1).value() / roomTLY) * Math.sqrt(Math.pow(theDPS.point(i).coordinate(1).errorMinus() / theDPS.point(i).coordinate(1).value(), 2) + Math.pow(roomTLYErr / roomTLY, 2)));
			theDPS.point(i).coordinate(1).setErrorPlus(theDPS.point(i).coordinate(1).errorMinus());
		}
		return theDPS;
	}


	/**
	 *  Normalize the LY to 1 at RoomT
	 *
	 *@param  theDPS   the LY DPS
	 *@param  roomTLY  The Value of the RoomT LY
	 *@return          The normalized DPS
	 */
	public static IDataPointSet normalizeToRoomT(IDataPointSet theDPS, double roomTLY)
	{
		return normalizeToRoomT(theDPS, roomTLY, 0.0);
	}


	/**
	 *  Normalize the LY to 1 at RoomT
	 *
	 *@param  theDPS  the LY DPS
	 *@return         The normalized DPS
	 */
	public static IDataPointSet normalizeToRoomT(IDataPointSet theDPS)
	{
		return normalizeToRoomT(theDPS, theDPS.point(0).coordinate(1).value(), theDPS.point(0).coordinate(1).errorMinus());
	}


	/**
	 *  Tells if an integer is even or not
	 *
	 *@param  theValue  The integer to evaluate
	 *@return           True if the integer is even
	 */
	public static boolean isEven(int theValue)
	{
		int theInteger = theValue / 2;
		double theDouble = theValue / 2.0;
		double theIntDoub = (double) theInteger;
		if (theDouble == theIntDoub)
			return true;
		else
			return false;
	}


	/**
	 *  Calculate the coincidence probability depending on the decay time values
	 *  and their fraction of light
	 *
	 *@param  amp        The array containing the amplitude value of each component
	 *@param  tau_us     The array containing the decay time value in us of each
	 *      component
	 *@param  numPhot    The number of photon on each channel
	 *@param  verbose    Display some information if true
	 *@param  tCoinc_us  Description of the Parameter
	 *@return            The coincidence probability
	 */
	public static double coincEfficiency(double[] amp, double[] tau_us, double tCoinc_us, double[] numPhot, boolean verbose)
	{
		if (verbose)
		{
			System.out.println("N0 = " + numPhot[0]);
			System.out.println("N1 = " + numPhot[1]);
		}
		double Ntr = 0;
		double Nr[] = new double[tau_us.length];
		for (int i = 0; i < tau_us.length; i++)
			Ntr += amp[i] * tau_us[i];

		for (int i = 0; i < tau_us.length; i++)
			Nr[i] = amp[i] * tau_us[i] / Ntr;

		double eff = 0;

		for (int i = 0; i < tau_us.length; i++)
			for (int j = 0; j < tau_us.length; j++)
				eff += Nr[i] * Nr[j] * tau_us[j] / (tau_us[i] + tau_us[j]) * Math.exp(-1.0 * tCoinc_us / tau_us[j]);

		double pCoinc = 1 - Math.pow(2 * eff, numPhot[0] * numPhot[1]);
		if (verbose)
			System.out.println("Pcoinc = " + pCoinc);

		return pCoinc;
	}


	/**
	 *  Calculate the coincidence probability depending on the decay time values
	 *  and their fraction of light
	 *
	 *@param  amp        The array containing the amplitude value of each component
	 *@param  tau_us     The array containing the decay time value in us of each
	 *      component
	 *@param  tCoinc_us  The length of the the coincidence window in us
	 *@param  numPhot    The number of photon on each channel
	 *@return            The coincidence probability
	 */
	public static double coincEfficiency(double[] amp, double[] tau_us, double tCoinc_us, double[] numPhot)
	{
		return coincEfficiency(amp, tau_us, tCoinc_us, numPhot, false);
	}


	///**
	// *  Removes the DataPoint with NaN on Y in an IDataPointSet
	// *
	// *@param  theDPS  The IDataPointSet with NaN
	// */
	//public static void removeNaN(IDataPointSet theDPS)
	//{
	//	IDataPointSet modDPS = af.createDataPointSetFactory(af.createTreeFactory().create()).create(theDPS.title(), 2);
	//	int j = 0;
	//	for (int i = 0; i < theDPS.size(); i++)
	//		if (theDPS.point(i).coordinate(1).value() != Double.NaN)
	//		{
	//			modDPS.addPoint();
	//			modDPS.point(j).coordinate(0).setValue(theDPS.point(i).coordinate(0).value());
	//			modDPS.point(j).coordinate(0).setErrorMinus(theDPS.point(i).coordinate(0).errorMinus());
	//			modDPS.point(j).coordinate(0).setErrorPlus(theDPS.point(i).coordinate(0).errorPlus());
	//			modDPS.point(j).coordinate(1).setValue(theDPS.point(i).coordinate(1).value());
	//			modDPS.point(j).coordinate(1).setErrorMinus(theDPS.point(i).coordinate(1).errorMinus());
	//			modDPS.point(j).coordinate(1).setErrorPlus(theDPS.point(i).coordinate(1).errorPlus());
	//			j++;
	//		}
	//
	//	theDPS = modDPS;
	//}

	/**
	 *  Shifts the X values of a P1D by a given value
	 *
	 *@param  p1d    The P1D to be shifted
	 *@param  shift  The shift value
	 *@return        A new P1D with the sa
	 */
	public static P1D shiftX(P1D p1d, double shift)
	{
		return shift(p1d, shift, 0);
	}


	/**
	 *  Shifts the Y values of a P1D by a given value
	 *
	 *@param  p1d    The P1D to be shifted
	 *@param  shift  The shift value
	 *@return        A new P1D with the sa
	 */
	public static P1D shiftY(P1D p1d, double shift)
	{
		return shift(p1d, shift, 1);
	}


	private static P1D shift(P1D p1d, double shift, int axis)
	{
		P1D shiftP1D = new P1D(p1d.getTitle(), false, p1d);
		double[] p1dX = shiftP1D.getArrayX();
		double[] p1dXL = shiftP1D.getArrayXleft();
		double[] p1dXR = shiftP1D.getArrayXright();
		double[] p1dXLS = shiftP1D.getArrayXleftSys();
		double[] p1dXRS = shiftP1D.getArrayXrightSys();
		double[] p1dY = shiftP1D.getArrayY();
		double[] p1dYU = shiftP1D.getArrayYupper();
		double[] p1dYL = shiftP1D.getArrayYlower();
		double[] p1dYUS = shiftP1D.getArrayYupperSys();
		double[] p1dYLS = shiftP1D.getArrayYlowerSys();

		for (int i = 0; i < p1dX.length; i++)
		{
			if (axis == 0)
				p1dX[i] += shift;
			else if (axis == 1)
				p1dY[i] += shift;
		}
		shiftP1D.clear();
		shiftP1D.fill(p1dX, p1dY, p1dXL, p1dXR, p1dYU, p1dYL, p1dXLS, p1dXRS, p1dYUS, p1dYLS);
		return shiftP1D;
	}

	///**
	// *  Shifts the X values of a IDataPointSet by a given value
	// *
	// *@param  dps    The IDataPointSet to be shifted
	// *@param  shift  The shift value
	// *@return        A new shifted IDataPointSet
	// */
	//public static IDataPointSet shiftX(IDataPointSet dps, double shift)
	//{
	//	for (int i = 0; i < dps.size(); i++)
	//	{
	//		shift(dps.point(i), 0, shift);
	//	}
	//}
	//
	//
	///**
	// *  Shifts the Y values of a IDataPointSet by a given value
	// *
	// *@param  dps    The IDataPointSet to be shifted
	// *@param  shift  The shift value
	// *@return        A new shifted IDataPointSet
	// */
	//public static IDataPointSet shiftY(IDataPointSet dps, double shift)
	//{
	//	for (int i = 0; i < dps.size(); i++)
	//	{
	//		shift(dps.point(i), 1, shift);
	//	}
	//}
	//
	//
	//private static void shift(IDataPoint pt, int axis, double shift)
	//{
	//	IMeasurement meas = pt.coordinate(axis);
	//	double val = meas.value();
	//	meas.setValue(val - shift);
	//}
}

